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I. INTRODUCTION 



There is a considerable freedom in formulating QCD on the lattice. This freedom is re- 
flected in the large number of actions tested and used in the quenched approximation. There 
are no miracles: good scaling, good chiral properties, theoretical safety and the expenses 
are in balance. Short of an algorithmic breakthrough one can expect to see in the future a 
plethora of full QCD simulations and results obtained with different formulations adapted 
to the physical problem as it happened in the quenched approximation. 

In this paper we discuss a full QCD algorithm for 2+1 light flavors with the parametrized 
fixed point action jl|, 3 [1Q| The lightest quark mass m ud which can be simulated is set only 
by the small chiral symmetry breaking caused by the parametrization error, the expenses of 
a full updating sweep are practically independent of m u d. 

Our updating procedure has no special problems in connecting different topological sec- 
tors nor in suppressing the topological susceptibility. The algorithm is exact and the action 
certainly describes QCD in the continuum limit. It is a partially global update where the 
pieces of the determinant are switched on gradually in the order of their expenses. The 
partially global update implies that this is a volume-squared algorithm which constraints 
the size of lattices one can cope with in the simulations. Actually, the expense of a specific 
section of the algorithm, the eigenvalue solver, increases even faster with the volume, but in 
the present simulations this part is not dominating. 

In the ongoing runs the target spatial sizes are L s « 1.2 fm and 1.8 fm with a resolution 
a ~ 0.15 fm. The target pion mass is below 300 MeV. On the 1.2 fm lattice we also ran with 
the smallest m u d quark mass which can be treated with our Dirac operator, basically with 
massless quarks. 

Actions which approximate the fixed point of a renormalization group transformation have 
been tested in detail in d = 2 and d — 4. The present form of the QCD fixed point action 
and the related codes are the result of a long development to which many of our colleagues 
contributed (see Ref . P| and references therein) . The algorithm discussed here is optimized 
and running on three different platforms (IBM SP4, PC Cluster, Hitachi SR8000). This 
paper is organized as follows. In Sect|n]we briefly describe the fixed point action. Sect lHIl 
summarizes the Partial Global Stochastic Update and its improvements in a general form 
while Sect II VI describes how the update is implemented in our simulation. In SectEl we 
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present some numerical results that illustrate the updating algorithm. 



II. THE ACTION 



A. The parametrized fixed point action 



The special QCD action which is the fixed point of a renormalization group transformation 
las several desirable properties It is a local solution of the Ginsparg- Wilson equation 

75 D + D 75 = D l5 2RD, (1) 



with a non-trivial local matrix R which lives on the hypercube |8j|. The quark mass is 
introduced as 

D{m) = D + mQ R - l -Dy (2) 

EqJT] guarantees that the Dirac operator is chirally invariant. Since it is the fixed point of 
a renormalization group transformation, it has no cut-off effects in the classical limit. The 
parametrized version of this action is an approximation which has been carefully tested in 
the classical limit and in quenched simulations 

The parametrized fixed point gauge action S g (U) |2| is a function of plaquette traces 
built from the original gauge links U, and from smeared links. The smeared link contains 
staples in an asymmetric way: the weights of staples which lie in, or orthogonal to the plane 
of the plaquette are different. The gauge action is a polynomial of smeared and unsmeared 
plaquette traces. The 5 non-linear and 14 linear parameters are fitted to the fixed point 
action. 

The Dirac operator |l|, Q| is constructed on smeared gauge configurations V{U). This 
smearing is local and contains links projected to SU(3). It is constructed using renormal- 
ization group considerations Il3j| and it reflects the discontinuous character of the chiral 

Q 

Dirac operator when the topological charge changes (3j. The Dirac operator has fermion 
offsets on the hypercube only. In Dirac space all the elements of the Clifford algebra enter. 
The structure of these terms is restricted by the symmetries C, P, 75-hermiticity, and cubic 
symmetry. The 82 free coefficients of this Dirac operator are determined by a fit to the fixed 
point Dirac operator. We note that also for the Chirally Improved Dirac operator - which 
has a similarly complex structure - a global update has been suggested in Ref. 14] 



B. The 2+1 flavor action 

Our goal is to simulate a 2 + 1 flavor system with quark masses close to their physical 
values. As usual we integrate out the the fermionic fields and write the action as 

S = j3S g (U) + uD(m ud )u + dD(m u d)d + sD(m s )s (3) 
~ /3S g (U) - \ndet 2 D(m ud ) - lndetD(m 5 ). (4) 

The 75-hermiticity of the Dirac operator ensures that det D(m) is real, detD(m) = 
detD^(m). If D were an exact solution of EqJT] then det D(m) would be positive for any 
m > 0. Due to parametrization errors our Dirac operator has no exact chiral symmetry 
and the requirement of det D(m) > puts a constraint on the quark mass m. For now we 
assume that this condition is satisfied. 

One can rewrite the action in an explicitly Hermitian form 

S = pS g (U) - In det (D\m ud )D(m ud )) - In det (j D^{m s )^ D{rn s )\ ■ (5) 

We want to emphasize that the appearance of the square root in EqlHl does not mean that 
we simulate a possibly non-local action. The action we simulate is given by the manifestly 
local form in EqOl 



III. THE STOCHASTIC UPDATE AND ITS IMPROVEMENTS 

The parametrized fixed point action contains many gauge paths and SU(3) projections. 
It is too complicated, if not impossible, to simulate with algorithms that would require the 
derivative of the action with respect of the gauge fields. For that reason we have adapted 
an updating method that requires only a stochastic estimate of the action at each step. 
(Partial) Global Stochastic Update was developed in lljj, based on an old suggestion 

nnn 

to simulate smeared link staggered actions. It was developed further in [lalllj 20]. 
Stochastic o^o^p dat,„ g algorithm have been used in different context by m any other 
groups 




I22J, [23j, m ake over the main points of the Global Stochastic Update 

but add several improvements to create an efficient updating algorithm. In the next section 
we will summarize the general ideas of the update, then discuss the specific improvements 
we have implemented. 
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A. The Partial-Global Stochastic Update 



In order to simplify the notation in this section we consider an action with the generic 
form 

S = pS g {U)-]ndetA*A. (6) 

Here A 1 A describes 1 or 2 flavors of massive fermions as discussed in Sect III Bl A Global 
Update proceeds in two steps: 

A: Update (a part of) the configuration U — > U' with the gauge action (3S g using 

Metropolis, over-relaxation or other updates. 

One could accept or reject (A/R) the proposed U' configuration with probability 

. / detA'tA'\ 

This procedure clearly satisfies the detailed balance condition but requires the evaluation 
of the fermionic determinant. This lengthy calculation can be replaced by a stochastic 
estimator as follows. We write the determinant ratio as a stochastic integral 

detA'U' , i 
where we introduced the notation 

Q = A'- 1 A. (9) 

B: For the stochastic accept /reject step we first create a Gaussian random vector 

i] with P(rj) oc exp(— ryrf). Now the new configuration is accepted with the 
probability 

P a s c t c och = min(l,e- A5 0, (10) 

where 

AS f = ^(fitfi - 1)77. (11) 

EqlTTI defines the stochastic estimator AS 1 /, the change of the fermionic action 
with fixed rj [l""" 
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The stochastic update satisfies the detailed balance condition 12, 3, and repeating steps 
A-B creates a sequence of gauge configurations with the proper probability distribution. 
For the 2+1 flavor system of EqJHlthe stochastic estimator contains two terms 

as/ = vU^ud - i)vu + JMn. - i)vs (12) 



with Q ud = D' ~ 1 {m ud )D{m ud ) and Q s = ^ D'^irris)^ D{m s ). 

The main ideas of the stochastic Monte Carlo update has been known for a long time but 
it has not been used in numerical simulations because of technical difficulties. In its original 
form the acceptance rate in EqUH] is infinitesimally small unless the configurations U and 
U' are nearly identical. The root of this problem is two-fold. 

On one hand, if the typical values of | logdet(D'/D)| are significantly larger than 1 the 
acceptance rate will be very small, even if the determinant ratio was calculated determin- 
istically. On the other hand, an additional suppression of the acceptance rate occurs due 
to the stochastic evaluation in the A/R step. To illustrate this consider the case when 
detD'/ det D = 1. Take a simple model for this situation with a Gaussian random variable 
P(x) oc exp(— (x — Xq) 2 /2a 2 ). The relation (e~ x ) = 1 implies x = <r 2 /2 hence for large a 
the acceptance rate is extremely small, ~ e~°" 2//2 <C 1. 

Note that the standard deviation of exp(— rf(£vQ — \)rj) is infinite if any of the eigenvalues 
of VL = A'~ 1 A is smaller than 1/2 However, the extremely small acceptance rate occurs 
much before this bound is reached. In the following we discuss four improvement steps which 
are essential to get an algorithm with a good acceptance rate. 



B. Improvements of the Global Stochastic Update 

In this section we discuss the different improvements we implemented to increase the 
effectiveness of the stochastic updating. In the reduction technique the UV part of the de- 
terminant is separated, its value is calculated non-stochastically and taken into account more 
frequently by intermediate A/R steps. This increases the acceptance rate in the stochastic 
A/R step significantly by reducing both problems mentioned above. The subtraction tech- 
nique separates the IR modes by calculating the eigenvalues and eigenvectors of the first few 
low eigenmodes of the Dirac operator. It acts analogously for the IR modes as the reduction 
for the UV part. The last two techniques, the relative gauge fixing and the determinant 
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breakup, are applied in the last, stochastic A/R step, and are aimed at reducing the fluctu- 
ations. The relative gauge fixing brings the configuration U' as close to U as possible. The 
determinant breakup rewrites the Dirac operator as product of operators. The stochastic 
estimator becomes a sum of independent terms and its fluctuation is reduced. 

The Global Stochastic Monte Carlo Update would not be effective without these im- 
provements. On large lattices even with the improvements one can update only a part of 
the configuration before evaluating the stochastic estimator making the algorithm to scale 
with the square of the volume. Nevertheless on moderate volumes we found the algorithm 
effective, allowing the dynamical simulation of light, even nearly massless, quarks with an 
action where the chiral breaking and lattice artifacts are small. 



1. The Reduction: 

The stochastic change of the fermionic action can be written as ASf = J^ii^i ~ l)ViVi 
where Ui are the real eigenvalues of the operator fl^fl and rjjrji = 0(1). While the eigenvalues 
of the Dirac operator are restricted to a compact region, LOi can vary between ~ m to ~ 1/m, 
though most of the eigenvalues correspond the the UV modes are 0(1). These UV modes 
contribute little to ASf individually, but there are so many of them that they dominate the 
fluctuations. To reduce the fluctuations we transform the Dirac operator D — > D r such that 
the UV modes of D r are condensed and thus the corresponding eigenmodes of fl are closer 
to unity. We choose D r such that the change in the determinant det(D/D r ) is calculable 
analytically (non-stochastically) . 

The effect of the reduction on the eigenvalue spectrum of D is illustrated in Figure [TJ The 
spectrum was calculated on a single 4 4 quenched gauge configuration U with a ~ 0.15 fm in 
|s|. The larger "Batman" like structure corresponds to the spectrum of the original Dirac 
operator. The "Batman ears" are mainly the consequence of the non-trivial R operator 
and are strongly reduced if one considers the spectrum of DR which is close to a circle. 
An overwhelming part of the eigenvalues are around the UV point (s, 0) in the complex 
plane where s ~ 2.8. The points to the right of the long dotted line represent 95% of the 
eigenvalues. With the reduction we attempt to move the eigenvalues of D/s close to 1. 

nnn 

Following Refs. |15j, [25j, |2j we define a reduced Dirac operator 



D^fle-^^K (13) 
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Figure 1: The eigenvalue spectrum of the Dirac operator on a single 4 4 gauge configuration. The 
points of the larger "Batman" like figure correspond to the original Dirac operator while the points 
of the smaller "wing-shape" structure in the center represent the corresponding eigenvalues of the 
reduced Dirac operator D r /s. Sections marked by A, B and C on the original spectrum are mapped 
to sections a, b and c after reduction. 

Choosing the coefficients q = (— the reduced operator is D r /s = l + 0((D/s — 1)" +1 ). 
The ratio of the determinants D r and D can be expressed in terms of traces of D 

detD r = detDe-^i c > Tr{D/s - lY 

= detDe-^i=i^ TrDi , (14) 

and can be calculated non-stochastically by evaluating Tr D k , k = l,n. The computing time 
and the complexity of the code to do the trace calculations increase rapidly as n increases. 
With our Dirac operator we decided to stop at n = 4 in the reduction. Working with D r 
is not much different from D. Both in multiplication and inversion the exponential term in 
EqUnican be approximated by a relatively low order polynomial. 

In Figure [U the smaller wing-shape object in the center corresponds to the eigenvalues 
of D r /s on the same gauge configuration as before. Sections marked as A, B and C on 
the original spectrum are mapped to sections a, b and c after reduction, illustrating how 
the eigenvalues are transformed. The origin is a stationary point. The main effect of 
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the reduction is condensing the UV modes. The points to the left of the short dotted 
line again represent 95% of the eigenvalues. This region is not very obvious in the figure 
since all its points are contained in a circle around (1, 0) with radius of about 0.01 (with 
the overwhelming part of the eigenvalues being much closer), showing the strength of the 
reduction. Accordingly, the UV fluctuations, that is the contribution of the UV modes to 
the stochastic estimator AS"/, is greatly reduced. 
At this point the action EqlHlcan be written as 

S = (3S g (U) + S uv 

- \ndetDl(m ud )D r (m ud ) - lndet D r (m s ) (15) 

where 

n n 

S uv = -2 £ a< Tr D\m ud ) - ]T on Tr D\m s ) (16) 
i=i i=i 

carries most of the UV part of the determinant. 

2. The Subtraction 

The reduction of the Dirac operator as discussed in the previous section effects mainly 
the non-physical UV part of the determinant. The subtraction that we introduce here 
deals with the low lying IR eigenvalues of the Dirac operator. The small eigenvalues of 
D' can create large u>i eigenvalues of Q = D' ~ X D. Besides suppressing such configurations 
in the (full QCD) equilibrium configurations, their presence produces large fluctuations in 
the stochastic estimator, reducing therefore the acceptance rate in the stochastic A/R step. 
By calculating some low lying eigenvalues (and the corresponding eigenvectors) one can 
take into account their contribution more frequently and deterministically, so they do not 
participate in the stochastic A/R step. 

Denote the right and left eigenvectors of the Dirac operator by 

Dv\ = \v\ , w[D = \w{ . (17) 
The eigenvectors w x can be chosen to fulfill the normalization condition 

= 5 \x' ■ (18) 
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In terms of the (non-hermitian) projector operators 



Px = v x w{ (19) 

which, due to EqUHl satisfy the relation P% — P\, the Dirac operator can be written as 

£ = £AP A . (20) 

The subtracted Dirac operator is defined by replacing a set of the lowest eigenvalues by the 
constant s 

D s = D + ^(s-\)P x . (21) 

low 

We assume that we subtract the complex conjugate pairs together. For an arbitrary analytic 
function f(D) the subtraction gives 

f(D s ) = /(£>) + £(/(*) - /(A))JV (22) 

low 

Subtracting the reduced Dirac operator of Eq. fT3l gives then 

D rs = De~^M D /^y + J2 (s - Ae-^^W-i)*) Px . 

low 

The small eigenvalues of the subtracted, reduced D rs operator are replaced by s while its 
UV part is condensed near s. The stochastic estimator of D rs has reduced fluctuations and 
reduced absolute value as well. 

The ratio of the determinants of D rs and D r can be calculated analytically using the 
relation 

det D r = det D rs JJ -e~ ^(V--i)' . (23) 

low ^ 



The inversion of D using that of D s is given by 



D; 1 ■ (24) 



The smallest eigenvalues of D are replaced by the constant s in D s therefore the conjugate 
gradient method converges faster for D s . 

At this point the action of Eqsl5l and H~5l can be written as 

S = (3S g {U) + S uv + S IR 

- In det Dl s (m ud )D rs (m ud ) - In det D rs (m s ) (25) 
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where 



S, 



IR 




low 



5. T/ie Relative Gauge Fixing 

Even if C/ 7 is a gauge transform of U, hence the determinant ratio is exactly 1, the 
eigenvalues of fl = D{U')~ l D{U) are in general different from 1 (only their product is 1). 
As a consequence the stochastic estimator in the A/R step can have large fluctuations, 
greatly reducing the acceptance rate 2^|. These fluctuations can be reduced significantly by 
gauge transforming U' so as to maximize Y,x,n ReTi^U^U^), i.e. by bringing U' as close 
to U as possible. (We have also tried to fix the gauge in both U and U' by some given 
"absolute" gauge fixing condition, but the method discussed above was more efficient.) 



4- The Determinant Breakup 

The reduction, subtraction, and relative gauge fixing result in significant improvement 
of the stochastic estimator. Further improvements can be achieved by writing the Dirac 
operator as the product of / terms 

A = A 1 x A 2 x .... x Ai . (27) 

The corresponding stochastic estimator is the sum of I terms 

= £771(^-1)77, (28) 
i=i 

with I stochastic rj vectors and = A\~ x Ai. If the eigenspectrum of the individual Ai 
operators is closer to a constant the fluctuation of the stochastic estimator is reduced. In 
Refs. (3, 0|- following a suggestion in Ref. the terms in EqEH were chosen to be 

identical, Ai = A 1 / 1 . While this choice does reduce the fluctuations, it creates / equally 
singular terms and requires the calculation of the Zth root for each of them. We found it is 

n 

more effective to generalize the mass shifting method of [25] and write the Dirac operator 
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as 



a^) AM " A(^) XA( ^' (29) 
where A(/i) = 7^ (A + fi) and /io = 0. The mass shift values /ij are chosen such that each 
term in ASf contributes approximately equally. 

The first term in EqEHis the easiest to analyze. At lowest order in Hi 

,f(njn, - Dm = -4 + % - % - v + <W& • (so) 

Unless the operators A and A' are close to each other, n%/A has to be small to control the 
stochastic fluctuations. Consequently \i\ has to be much smaller than the smallest eigenvalue 
of A. Later terms allow larger change in the shift masses fii- The last mass of the series, 
Hi, is chosen such that the stochastic estimator of the single operator A(fii) is comparable 
to the previous terms. It is interesting to note that the leading term in EqEDl vanishes for 
a Ginsparg- Wilson Dirac operator with R = const. However it is not zero in our case. 

In practice we combine the reduction, subtraction, relative gauge fixing and the determi- 
nant breakup. Only the last two terms of EqEHlare treated stochastically. For the degenerate 
u and d quarks we can have 

A = D rs (m ud ) , 
A(fi) = ( £l^+£ ) , (31) 

where the subtraction is defined as in Eqj22j Each term in the stochastic estimator requires 
two multiplications and two inversions by A(/i). For the inversion a standard conjugate 
gradient or its variant can be used. 

For the s quark the situation is slightly more complicated. The A operator contains a 
square root operation 



A = \/D rs (m s ) , 



A(h) 



( D r (m s ) + ji 



(32) 



Again, each term in the stochastic estimator requires two multiplications and two inversions 
by A(/i). For both of these we approximate the square root operator by a polynomial 
series. Polynomials have been used to approximate both positive and negative roots of 
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Dirac operators before 1 10l. I27L 29]. Our situation is different because the operator D r (m s ) 
is complex. The case of optimal polynomials for a complex spectrum has been studied in 
Fortunately the strange quark mass is sufficiently heavy and a Taylor expansion in 
(D(m)/s — 1) works well. 



IV. THE ALGORITHM WITH 2+1 FLAVORS 

We describe now the algorithm which has been coded, tested and optimized for different 
platforms. The algorithm starts with a partially global gauge update which is followed by 
several accept/reject steps, where parts of the determinant are switched on gradually in the 
order of their expenses. It is convenient to rewrite the action of Eql2H]in a different form 

S = ((3 + 5(3)S a (U) 
+ [S 9 UV -5(3S 9 {U)\ 

+ [Suv - S 9 UV + Sj^ T ] (33) 
+ [S IR - S^ pT - In det Dl s (m ud )D sr (m ud ) - In det D rs (m s )]. 

The meaning of the different terms will be explained in the rest of this section. 



A. Gauge update 

The gauge update is a standard Metropolis/over-relaxation local update with the fixed 
point gauge action at coupling (3 e R = (3 + 5(3. Here 5(3 (added at this point and subtracted 
later) approximates the shift of the gauge coupling due to the determinant and helps to 
generate configurations with lattice spacing a that is close to the target value already at 
this step. 4n p gauge links, originating from n p consecutive lattice sites, are updated with 
Metropolis and then the same links with over-relaxation in a reversible sequence. This 
combination of updates we shall call a 'double update' in the following. 

In the first test runs our target lattice spacing is a ~ 0.15 fm and n p is 128 and 144 on 
the 8 3 x 24 and 12 3 x 24 lattices, respectively. In order to get the resolution a close to the 
target value we had to tune the coupling (3 repeatedly. The figures in this paper refer to the 
choice (3 = 3.15. 
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Figure 2: The correlation between ASjjv and A5 9 . The slope predicts the optimal value of 5(5 « 



B. The 1st accept/reject step 

The gauge configuration created as discussed above is accepted/rejected (A/R) with 
the action S uv — 8j3S g (U), where S uv is a good gauge approximation to the reduction 
contribution Suv of EqUHl The function S uv is represented by different gauge loops with 
fitted coefficients on the smeared configuration V(U). This smearing is the same which was 
used in the parametrized Dirac operator (Sect lH A"} . Calculating S uv is fast and can be 
done without building up the Dirac operator. The deviation between S uv and the exact 
reduction Suv will be corrected in the 2nd accept/reject step below. The parameter 6(3 is 
chosen to maximize the acceptance rate in this step. Figure [2] shows the correlation between 
AS g , the change of the gauge action, and ASw, the change of the contribution from the 
reduction, for a set of configuration pairs {U,U'}. From the slope 6(3 = —0.15 seems to 
be a reasonable choice. It is interesting to note that for our action the introduction of 
the determinant increases the gauge coupling. In the usual Wilson and staggered fermion 
simulations this shift is larger and in the opposite direction. Since the reduction contribution 
is large for distant configurations and the —6(3S g term cancels it only approximately, we have 
to keep the number of updated links 4n p in the gauge update modest in order to get a good 
acceptance rate in this 1st accept /reject. The combination of steps in Sects IIV Al and HV Bl 
is repeated N\ times. 



0.15. 
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Figure 3: The approximation of ASw in terms of gauge loops, AS UV . This error is corrected in 
the 2nd A/R step. 

With the n p value quoted before, the 1st acceptance rate is above 0.5. In the running 
simulations N\ = 28, i.e. 4 X n X N\ ~ 15k links are double updated before the 2nd 
accept/reject step. 

C. The 2nd accept /reject step 

The cycle of repeated steps in Sects ITVAl and ITV Bl is followed by a 2nd accept/reject 
decision. In this step the Dirac operator is built on the U' competitor configuration, the 
traces are calculated for the exact reduction, and a certain number of the lowest eigenvalues 
and eigenvectors are determined. The proposed configuration is accepted/rejected with the 
action (Suv ~ S uv ) + S^ pr . The first term corrects the small error we made in the 1st A/R 
step in approximating the traces in Suv with gauge loops. This error is typically small as 
shown in Figure El where the change AS UV calculated in the gauge approximation is shown 
as the function of its exact value. (The action differences plotted in this figure are taken 
between configurations which are offered to the 3rd A/R step, as a result of several 2nd A/R 
steps. These large values of O(10) would cause a very small acceptance rate in the 3rd step, 
had we not taken into account this contribution more frequently in the 2nd step.) 

The last term Sj R pv is an approximation to the contribution of the low lying eigenvalues to 
the determinant Sir in EqEBl ^ev eigenvalues and the corresponding Sir are calculated for 
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m U( i using an Arnoldi eigenvalue finder. The eigenvalues for m s are determined from these, 
using leading order perturbation theory. (Due to the presence of R in the Ginsparg- Wilson 
relation, EqJH the reaction of the eigenvalues on changing the quark mass is not a simple 
shift.) This approximation is very good, the error is typically O(10~ 4 ) . The combination 
of the steps in Sects lIV Al IIV Bl and IIV Cl is repeated N 2 times. 

At present we calculate n ev = 48 eigenvalues, which include all the eigenvalues with ab- 
solute value below « 0.40 and « 0.19 on the 8 3 x 24 and 12 3 x 24 lattices, respectively. 
The reduction shifts these values further to around ~ 1.0 and ~ 0.5. We use an Arnoldi 
routine from the publicly available PARPACK package. The application is far from optimal 
in our case. Even though the eigenvalues and eigenvectors are calculated on very similar 
configurations and change little from step to step, the routine cannot use this fact and cal- 
culates each eigenvalue set basically independently. The internal calculations of the package 
are also expensive, the multiplication of a vector by the Dirac operator does not dominate 
it. These problems limit the number of times the 2nd A/R step is repeated. We can afford 
N 2 = 6 repetitions and the acceptance rate of the 2nd A/R step is around 0.65. Overall 
about ~ 90k links are double updated before the 3rd A/R step. 

D. The 3rd accept/reject step 

The cycle described above is followed by a final, stochastic accept/reject step with the 
action Sir — 5"^ pr — In det Dl s (m u d)D rs (m u d) — lndet D rs (m s ). The first part corrects the 
small error we made in calculating the contribution of the low lying eigenvalues of D(m s ) 
to the determinant in the 2nd A/R step. For that we determine the low lying spectrum of 
D{m s ) on the competitor configuration U' which is first relative-gauge-fixed with respect to 
U. The second term gives the stochastic estimator of the subtracted, reduced, 2+1 flavor 
determinant. For the light quarks we break up the determinant into 76 terms. The lowest 
eigenvalue of the reduced subtracted Dirac operator, D rs /s is around one and the first few 
mass shifts have to be much smaller than that (see Sect JIII B~4jl . We chose A/i, = yu i+1 — /ij = 
0.01 for i = 1 — 20. The later A/x values are considerably larger. Due to the subtraction of 
the low lying eigenmodes the conjugate gradient iteration converges relatively fast, in about 
70 steps at the lowest mass shifts and in 5-10 steps at the largest, each step requiring two 
Dxv Dirac operator multiplications. The exponential term for the reduction requires about 
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20 D x v multiplications. For m s the determinant breakup has 38 terms and the smallest 
shift is A/i = 0.02. The square root and its inverse of the reduced, subtracted Dirac operator 
is approximated by their Taylor series in (D/s — 1). In the smaller mass shift region we use 
250-300 order polynomials, for the larger mass shift values this reduces to order 30-40. We 
expect that this section of the code could be significantly improved. 

With this closing accept/reject step the algorithm becomes exact. The steps in Sects. 
IIV A[ IIV B\ IIV C[ and IIV Dl are repeated and the accepted configurations that went through 
all three filters form a Markov chain corresponding the parametrized fixed point action. The 
acceptance rate of the 3rd A/R step is approximately 0.4. 

E. Optimization and performance 

The code is optimized on three different platforms (IBM SP4, PC cluster and Hitachi 
SR8000). The dominating numerical step is the multiplication of a vector by the Dirac 
operator, D x v, which can be effectively parallelized on all the platforms. The effectiveness 
of the internal manipulations of the PARPACK package, however, is very sensitive to the 
architecture. It would be very preferable to replace this part of the code by a QCD specialized 
piece. 

The stochastic estimator (in the 3rd accept /reject) requires ~ 21k and ~ 26k D x v 
multiplications on the 8 3 x 24 and 12 3 x 24 lattices, respectively. To calculate the first 48 
eigenvalues/eigenvectors of the Dirac operator requires ~ lk and ~ 2k D x v multiplications 
on the smaller and larger lattices, respectively. 

At certain stages of the calculation the processors are divided in two groups and work on 
the gauge and Dirac part of the code independently. They are joined, however, to calculate 
the stochastic estimator together. 

V. PRELIMINARY RESULTS 

In our first set of test runs, starting from two different a ~ 0.15 fm quenched configura- 
tions, we generated about 400 + 600 8 3 x 24 configurations, each separated by a full cycle of 
updates and A/R steps as described in Sects lIV AllV Dl We chose the run parameters, based 
on earlier quenched runs, as (3 = 3.15, 5/3 = —0.15, m u d = 0.017 and m s = 0.095. Figure 21 
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sweeps 



Figure 4: The equilibration of the plaquette as the function of updating sweeps in a 2 + 1 flavor 
dynamical simulation. The two series started from two different quenched gauge configurations. 
One sweep on the horizontal axis corresponds to a double update of the whole lattice. 

shows the equilibration of the plaquette as the function of updating steps for our two sets. 
The units on the horizontal axis are given in sweeps that correspond to a double update 
of the whole lattice which corresponds to about 5 full cycles. We determined the lattice 
spacing from the static potential as a = 0.14(l)fm, which gives the spatial size L s « l.lfm. 
On this rather small volume the hadron spectrum shows large finite volume effects. Instead 
of the pion mass, which is dominated by the volume, we estimate the quark mass from the 
eigenvalue spectrum of the Dirac operator. 

The left panel of FigureElshows the first 48 low energy eigenvalues on 50 a ~ 0.15 fm pure 
Yang-Mills configurations with m q = 0.017. This quark mass corresponds to m n ~ 300 MeV 
pions in the quenched approximation. The right panel shows the first 48 low lying eigenvalues 
on 50 equilibrated dynamical configurations. The eigenvalues of a massless chiral Dirac 
operator lie on circle if R = const. Our Dirac operator has a non-trivial R in which case 
one knows crude bounds only:the eigenvalues should lie between two circles touching each 
other at the origin. The full spectrum on a 4 configuration in Fig. Q] gives more information. 
A small quark mass shifts the eigenvalues to the right. The scattering of the eigenmodes 
characterize the chiral symmetry breaking of our approximate Dirac operator. 

While the scatter of the eigenmodes on the left panel is not negligible, its scale is small 
(compare to the whole spectrum of Fig. [T]). In the quenched hadron spectrum calculations 
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of Refs. |s, IS, 3] at similar parameters no exceptional configurations were observed which 
is consistent with the quenched spectrum in Fig. [21 here. 

There are several new features we can identify on the right panel that corresponds to the 
dynamical configurations. First we observe that the eigenmodes are shifted somewhat to 
the left. Their value suggests a nearly zero physical quark mass indicating a small additive 
mass renormalization of 8m ~ —0.015. On the quenched configurations 8m is practically 
zero. As we based our parameters on the quenched spectroscopy, by accident we simulated 
an approximately massless dynamical quark system. Since in the 2nd A/R step of the 
update (Sect. IIV CI) we subtract the low lying eigenmodes, this did not cause any increase 
in computing time. 

On its own a small mass renormalization is not a problem. It is the fluctuations of 
the eigenmodes beyond 8m that create exceptional configurations. These fluctuation are 
suppressed on the dynamical configurations as compared to the quenched case. The two 
panels of Figure [U contain the same number of eigenvalues and it is apparent that the Dirac 
operator on the dynamical configurations is much more chiral than on the quenched ones. 
This unexpected benefit is the effect of the fermionic determinant in the Boltzmann weight 
and shows that its presence enhances the configurations where our parametrization of the 
fixed point Dirac operator works better. 

The very small eigenmodes, those with |A| < 0.1, are completely missing from the right 
panel. This is the consequence of the suppression of the low eigenmodes by the determinant. 
The gauge update of Sect IIV Al creates configurations with real eigenmodes and some of 
these are accepted by the A/R steps. One of these modes is present on the right panel of 
Figure El at A « 0.3. However as these real eigenvalues move toward zero their determinants 
become small and the configurations are eventually replaced by configurations without real 
eigenmodes. On large volumes the small complex eigenmodes are not completely suppressed 
but on these small volumes even those are missing. 



VI. CONCLUSION 



In this paper we discussed the Partial Global Stochastic Update method and its use to 
simulate dynamical fixed point fermions. The original method is combined with several 
improvement techniques. By separating and controlling both the UV and IR modes of 
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Figure 5: The low lying eigenvalue spectrum of the Dirac operator on 50 pure Yang-Mills (left) and 
50 equilibrated dynamical (right) configurations. Both panels correspond to 8 3 x 24, a m 0.14fm 
lattices with the same lattice quark masses. 

the Dirac operator and the fluctuations of its stochastic estimator we found the update 
efficient on moderate volumes even near the chiral limit. To illustrate the algorithm we 
presented some preliminary results on 8 3 x 24, a ~ 0.14 fm lattices with 2+1 flavors with an 
approximately massless light doublet. 
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